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Dynamical behavior of steady granular flow is investigated numerically in the inelastic hard 



flow regimes, i.e., the collisional flow and the frictional flow. In the collisional flow, the hard sphere 



value and the total contact time fraction with other particles goes to zero. For the frictional flow, 
however, we demonstrate that the collision rate diverges as the power of the particle stiffness so 
that the time fraction of the multiple contacts remains finite even in the hard sphere limit although 
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sphere limit of the soft sphere model. We find distinctively different limiting behaviors for the two 
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limit is straightforward; the number of collisions per particle per unit time converges to a finite 



the contact time fraction for the binary collisions tends to zero. 



S 




'''Electronic address: naka4scp@mbox.nc.kyushu-u.ac.jp 



1 



I. INTRODUCTION 



The interactions between grains in flowing granular material are roughly classified into 
two types; the impulsive contact (collision]with the momentum exchange and the sustained 
contact with the transmission of forces The flow in which the impulsive contact is 
dominant is called collisional flow, while the flow where the sustained contact dominates is 
called frictional flow. These two types of flow may be found in a simple geometry such as 
granular flow on a slope. The grains stay at rest when the inclination angle is too small. If 
the inclination exceeds a critical angle, the material starts flowing frictionally at first. The 
flow becomes collisional when the inclination is large enough. 

As for the collisional flow of granular material, its dynamics has some analogy with 
molecular fluid, and the kinetic theories based on inelastic binary collisions of particles hold 

n 

to some extent [2J. On the other hand, the frictional flow is drastically different from the 
molecular fluid, and we have little understanding on it. Many models have been proposed 
for dense granular flows: For example, some models take into account the effect of non- 
local force transmission which comes from the network of contacting grains |3[. In the 
experiments, however, it is difficult to specify the nature of sustained contact in the dense 
flow. 

For the simulations of granular dynamics, the following two models have been commonly 
used, i.e. the inelastic hard sphere model and the soft sphere model. 

In the inelastic hard sphere model, the particles are rigid and the collisions are thus 
instantaneous, therefore its dynamics can be defined through a few parameters that charac- 
terize a binary collision because there are no many-body collisions ^J. The model is simple 
and there are very efficient algorithms to simulate it [oj], but it describes basically only the 
collisional flow || because the sustained contact is not allowed. It is also known that the 
system often encounters what is called the inelastic collapse J?l 0, 0, 13 ; the infinite number 
of collisions take place among a small number of particles in a finite time, thus the dynamics 
cannot be continued beyond that point without additional assumption. 

On the other hand, in the soft sphere model, that is sometimes called the discrete element 
method (DEM) in the granular community P, Q], the particles overlap during collision and 
the dynamics is defined through the forces acting on the colliding particles. Collision takes 
finite time, and not only binary collision but also many-body collision and sustained contact 
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between particles are possible, therefore, both the collisional and the frictional flows may 
realize in the model. Many researches have been done on granular flow down a slope using 
the DEM in both the collisional regime and the frictional regime [12I IkI Q| . 

In actual simulations, however, the stiffness constants used in the soft sphere model are 
usually much smaller than the one appropriate for real material such as steel or glass ball 
[15I ] because of numerical difficulty. Therefore, some part of sustained contact in simulations 
may be decomposed into binary collisions if stiffer particles are used. It is also not clear 
that the frictional force in the sustained contact may be described by the same forces with 
the one suitable for the collisional events. 

It is, therefore, important to examine how the system behavior may change as the stiffness 
constant increases in the soft sphere model, or in other words, how the soft sphere model 
converges to the inelastic hard sphere model in the infinite stiffness limit. In this paper, we 
present the results of numerical simulations of the granular flow using the soft sphere model, 
and investigate the system behavior when we change the stiffness constant systematically 
with keeping the resulting restitution constant unchanged. 

In Sec. m we briefly summarize an inelastic hard sphere model that is used for collisional 
granular flows. Then we introduce a simple soft sphere model for granular material and 
discuss how we take the hard sphere limit keeping the restitution coefficient constant. In 
Sec. IIII1 the simulation results are shown. At first, the stiffness dependence of the steady 
state of a single particle rolling down a slope is presented to see how the inelastic collapse 
appears in the hard sphere limit. The collisional flow and the frictional flow are examined. 
We find that the interactions between particles in the collisional flow smoothly converge to 
binary collisions of inelastic hard spheres, while those in frictional flow shows non-trivial 
behavior; the behavior is also different from the "inelastic collapse" in the single particle 
system. The summary and the discussion are given in Sec. IIVI 



II. HARD SPHERE MODEL AND SOFT SPHERE MODEL 

In this section, we introduce the hard sphere model and the soft sphere model for granular 
material. For simplicity, we consider a two dimensional system and grains are modeled by 
two dimensional disks. 

After showing the correspondence of parameters in the soft sphere model to those of the 
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hard sphere model, we briefly summarize the phenomenon called inelastic collapse, which is 
the singular behavior in the inelastic hard sphere system. At the last part of this section, 
we discuss the limiting behavior of the soft sphere model that corresponds to the inelastic 
collapse of the hard sphere model. 



A. Inelastic hard sphere model 

In hard sphere models, collisions between particles or between a particle and a wall are 
considered to happen instantaneously, and its dynamics is defined by the binary collision rule. 
We consider the collision rule in terms of the normal and tangential restitution coefficients 
and the sliding friction 0, Q] . 

Let us consider a collision between the two spheres i and j of the diameters Cj and err- 
and the masses mj and rrij at the contact positions r*j and rj, respectively. Prior to the 
collision, the disks have velocities c$ and Cj and angular velocities u>i and uij. Then the 
relative velocity of the point of contact Vij is given by 

vn = (cj - Cj) + n x + ^u-j , (1) 

where the normal vector n = r,ij/\rij\ = (n x , n y , 0) with = r, — rj. 

If v[j denotes the post-collisional relative velocity, the collision rule for normal direction 

is 

{V>'v'ij) = -e{n-Vij), (2) 

where e is the normal restitution coefficient with < e < 1. 

In the case without sliding, the collision rule in the tangential direction is given by 

(n x w y = -p{n x Vij ), (3) 

where (3 is the tangential restitution coefficient with — 1 < (3 < 1. The sliding is taken 
into account so that the tangential component of impulse does not exceed fi\n ■ J\ with the 
Coulomb friction coefficient /x, where J is the momentum change of the particle i through 
the collision. Namely, when the momentum change of the particle i through the collision 
rules (J2J) and (J3J) is J nos , then (n x v[-) is determined so that 

\n x J\ = min(/i|n • J|, \n x J nos \) (4) 



is satisfied. The detailed description of the rule is given in Refs. 



B. Soft sphere model: DEM 



The DEM, or the soft sphere model jj, is often used to simulate the dynamics of 



granular materials. In the present work, we adopt the two dimensional one with the linear 
elastic force and dissipation. When the two disks % and j at positions r, and Vj with velocities 
Cj and Cj and angular velocities u)i and ojj are in contact, the force acting on the particle i 
from the particle j is calculated as follows: The normal velocity v n , the tangential velocity 
Vt, and the tangential displacement u t are given by 

v n = n-Vij, v t = t-Vij, (5) 
u t = v t dt, (6) 

Jt 

where is given in Eq. ffT|) - and t = (—n y , n x , 0). Here, to is the time when the two particles 
start to be in contact. Then the normal force FJ] and the tangential force Ff, acting on the 
particle % from the particle j are given by 

FTj = 2Mk n (?i±?i _ | r . .\ j _ 2M Vn v n , (7) 

F\ 3 = min(|^|,/i|F n |)sign(^) (8) 

with 

h t = -2Mk t u t - 2Mrj t v t , (9) 

where k n and fct are the elastic constants, r] n and rjt are the damping parameters, [i is the 
Coulomb friction coefficient for sliding, and M = miirij/ '{mi + rrij) is the reduced mass. 



C. Hard sphere limit of soft sphere model 

When we adopt the linear dependence of the elastic and viscous force on the overlap and 
the normal relative velocity as in Eq. (JJJ), we can calculate the duration of contact r c and 
the restitution coefficient for a normal binary collision e ^?|; they are given by 

Tc= * 2 (10) 



and 



exp [ nVn ) , (11) 



respectively. 

Neglecting the sliding friction and the variation of n during the contact, the half period 
of the oscillation in the tangential displacement u t is estimated as 01 



7T 



(12) 



for the two dimensional disks with the moment of inertia Jj = mjCrf/8. Here we choose the 
parameters k t and rj t so that the relation r s = r c , or, 



y/6k t - 977 t 2 = ^/2k 



n 'In 



(13) 



is satisfied. Under this condition, the tangential restitution coefficient is given by 



/3 = exp(-3^r c ). 



(14) 



Equations (|10j) and (|11|) can be rewritten as 

1k n 



Vn 



(vr/lne) 2 + 1 



1/2 



7T 2 + (lne) 



1/2 



(15) 
(16) 



Using Eqs. ([13 )1 - ()lb |) . we can take the inelastic hard sphere limit of this model for given e 
and (3 by taking the k n — ► 00 limit; r] n , r] t , and k t are given by Eqs. (fI5j). JUJl, and 
and the duration time of collision r c goes to zero as Eq. fjl6j) . 



D. Inelastic collapse 

It is well known that the inelastic hard sphere system can undergo the inelastic collapse, 
i.e., the phenomenon where infinite collisions take place within a finite period of time 0. 
The simplest example is the vertical bouncing motion of a ball under gravity, but it has 
been shown that the inelastic collapse also occurs in higher dimensional systems without 
gravity QSSQ. 

The inelastic collapse never occurs in the soft sphere model because of the finite length 
of the collision time. However, it is worth to discuss what will happen in the soft sphere 
model in the simple situation where the inelastic collapse occurs in the hard sphere model. 
First, let us consider the vertical bouncing motion of a soft ball under gravity. We assume 
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the same force law between the ball and the floor as in Eqs. ((7J) and (JHJ), except that we 
replace 2M by the mass of the ball. While the ball and the floor are in contact, the equation 
of motion for the overlap of the ball and the floor, z, is given by 

z + k n z + rj n z = g, (17) 

with the acceleration of gravity g, where the dots represent the time derivatives. One can 
show that, if the impact velocity V\ is below a critical value v c , which is 0(1/ \fk~n) for large 
enough k n , the ball stays in contact with the floor; for Vi ^> v c , the restitution coefficient 
e w can be considered as a constant. Therefore, for a given initial impact velocity t>o, the 
number of necessary collisions n c for the ball to stay in contact with the floor is roughly 
estimated by the condition e™ c t>o ~ v c , namely, n c ~ In (v c /vq) / \ne w . Because v c ~ 1/v^n 
for large k n , n c behaves as 

n c oc In k n + const. (18) 

in the hard sphere limit. Thus n c diverges logarithmically when k n — > oo, which corresponds 
to the inelastic collapse due to gravity. 

In the case without gravity, the inelastic collapse results in the "many-body collision" in 
the soft sphere model. For example, we consider three soft spheres in one dimension and 
assume that the binary collision can be approximated by the collision law with a constant 
restitution coefficient. In the situation where a particle goes back and forth between the two 
particles approaching each other, the inelastic collapse may take place in the hard sphere 
model ; three balls lose relative velocity completely in the limit of infinite collisions. In the 
soft sphere model, however, when the interval between two collisions becomes smaller than 
the duration of contact r c , three balls are in contact at the same time, namely, the three 
body collision occurs; then they will fly apart. The number of collisions before the three 
body collision will also diverge logarithmically in the hard sphere limit because r c oc 1/ \fk^. 

On the other hand, one should note that a many-body collision in the soft sphere model 
does not necessarily result in the inelastic collapse in the hard sphere limit. Actually, in 
most of the cases, a many-body collision will be decomposed into a set of binary collisions 
in the hard sphere limit. 
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III. SIMULATION RESULTS 



In this section, we investigate the stiffness dependence of granular material on a slope 
in the following three situations: (i) a single particle rolling down a slope, (ii) the dilute 
collisional flow, and (iii) the dense frictional flow. We focus on the steady state in each 
situation and compare the simulation data with changing k n systematically. The particles 
are monodisperse in (ii), while they are polydisperse in (iii) in order to avoid crystallization. 

In the simulations, the parameters have been chosen to give e = (3 = 0.7, \x = 0.5 in the 
hard sphere limit. All values are non-dimensionalized by the length unit a, the mass unit m, 
and the time unit -J ojg. Here, a is the diameter of the largest particle in the system and m 
is the mass of that particle. The second order Adams-Bashforth method and the trapezoidal 
rule are used to integrate the equations for the velocity and the position, respectively [18]. 
Note that the time step for integration, dt, needs to be adjusted as r c becomes smaller. All 
the data presented in the paper are results with dt = min(r c /100, 10~ 4 ). We have confirmed 
that the results do not change for dt < r c /100 by calculating also with dt = r c /50 and 
dt = r c /200 in the case of the single particle. 



A. A single particle rolling down a bumpy slope 

Let us first consider a single particle rolling down a bumpy slope. It is known that the 

particle shows the steady motion for a certain range of the inclination angle 9 [19]. In the 

simulations, we make the boundary rough by attaching the same particles with the rolling 

one to the slope with the spacing 0.002a (see Fig. For the chosen parameters with 

the normal stiffness kn = 2~ 1 x 10 5 , the range of 9 for which steady motion is realized is 
\\ 

0.11 < sin# < 0.14 [13| . Here we fix the inclination angle to sin# = 0.13. Figure El shows 
the time evolution of the velocity in the y direction, v y , with k n = 2~ 1 x 10 5 (solid line) and 
k n = 2 17 x 10 5 (dashed line). It is shown that v y behaves periodically; this period (At ~ 4) 
corresponds to the period for the particle to get past one particle at the floor. The period 
hardly depends on the stiffness. 

In Fig. El (a), the k n dependence of the time averaged kinetic energy E is shown. It is 
shown that the energy is an increasing function of l/k n in the softer region (l/k n > 10~ 7 ), 
but no systematic l/k n dependence of E exists for l/k n < 10~ 7 . The average collision 
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rate (number of collisions per unit time) between the slope and the particle , N w , shows 
logarithmic dependence on l/k n in the whole region (Fig. 01(b)). From Fig. 0K C ) ; we also 
find that the average contact time fraction with the slope, t w , is a decreasing function of k n 
in the soft region, but it seems to approach a constant value for large enough k n . 

The logarithmic k n dependences of N w and the constant t w in large k n region agree with 
our previous analysis of "inelastic collapse under gravity" in the soft sphere model in Section 
III Dl The motion of the particle in one period is as follows; when the particle comes to a 
bump (a particle attached to the slope), the particle jumps up, bounces on the bump many 
times, loses the relative velocity, and finally rolls down keeping in contact with the bump. 
Therefore, the contact time fraction has a finite value even in the hard sphere limit due to 
the rolling motion at the last part. N w increases logarithmically in the hard sphere limit as 
has been expected from Eq. (JTSJl . 



B. Collisional flow 



Next we consider the steady state of the collisional flow. The system considered is shown 
in Fig. El The particles are monodisperse, and the slope is made rough as in the single 
particle case. The periodic boundary condition is adopted in the flow direction (x direction). 
The length of the slope is L = 50.1 and the number of the particles attached to the slope 
is 50. The number of flowing particles is also 50, namely, the number of the particles per 
unit length along the slope is about 1. The inclination angle is set to be sin# = 0.45. The 
initial configuration of particles is the row of 50 particles at rest with regular spacing in the 
x-direction, but each particle is at random height in the y-direction. After a short initial 
transient, if the total kinetic energy fluctuates around a certain value, we consider it as the 
steady flow. All the data were taken in this regime and averaged over the time period of 
1500. 

As can be seen in the snapshot, Fig. 01 the particles bounce and the number of particles 
in contact is very small. In Fig. 01 (a), the averaged kinetic energy per one particle, E, 
is shown. Though E becomes larger as the particles become softer in l/k n > 10 -5 , the 
systematic dependence of E on k n disappears for large enough k n (l/k n < 10~ 6 ) 2^|. The 
^/-dependence of the average number density and the flow speed in this region are also shown 
in Fig. 01(a) and (b). 
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In Fig. 03(b), the average collision rates between particles, N c (filled circles), and between 
particles and the slope, N w (open circles), per particle are plotted vs l/k n . Both of them 
stay roughly constant in the region where E is almost constant, l/k n < 10~ 6 . On the other 
hand, the average time fractions during which a particle is in contact with other particles, 
t c (filled circles), and in contact with the slope, t w (open circles), decrease systematically as 
k n increases as shown in Fig. 03(c). Both the solid and the dashed lines are proportional to 
namely, t c decreases in the same manner as r c in Eq. (|16p. Actually, the collision 
time fractions t c and t w converge to N c r c and N w t w (t w is the duration of a normal collision of 
a particle and the floor 2l| ) , respectively, in the hard sphere limit. The differences t c — N c r c 
and t w — N w t w are plotted in Fig. El (d) to show that they go to zero very rapidly upon 
increasing k n . This means that the interactions of the soft sphere model in the collisional 
flow regime converge to those of the inelastic hard sphere model with binary collisions. 

If we look carefully, however, in the large k n region in Fig. 0(c), we can see slight deviation 
of t w from the dashed lines; it decreases slower than 1 / \fk^. This tendency may indicate that 
t w remains finite in the hard sphere limit: Actually, in the event-driven simulation of the 
hard sphere model, we always found the inelastic collapse as long as the restitution constant 
between a particle and the floor is less than 1. This suggests that N w should diverge and t w 
should remain finite in the hard sphere limit because of the inelastic collapse due to gravity. 
The slight deviation of t w from the dashed line in Fig. 03 (c) may be a symptom of it, while 
we cannot see the logarithmic divergence in N w . 



C. Prictional flow 



Now we study the steady state of the dense frictional flow. We adopted the flat boundary 
with slope length L = 10.02 and imposed the periodic boundary condition in the flow 
direction. In order to avoid crystallization, we used the polydisperse particles with the 
uniform distribution of diameter from 0.8 to 1.0. The number of particles in the system is 
100. The inclination angle 9 is set to be sin# = 0.20 and the initial condition is given in 
the similar way with the collisional flow. The ten rows of ten particles with regular spacing 
in the x-direction are at rest with large enough spacings between rows in the y-direction 
so that particles do not overlap; only the particles in the top row are at random heights in 
y-direction. 
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It turns out that the steady state is not unique and fluctuation is large. Depending on 
the initial conditions, the particles may flow with a different value of average kinetic energy, 
or in some cases, the whole system stops. We presume that this is mainly because the 
system is not large enough: Grains form a layer structure as in Fig. and the flow velocity 
strongly depends on the configuration of particles in the bottom layer. The time evolutions 
of the kinetic energy per particle, E(t), of three samples with different configurations in the 
bottom layer are shown in Fig|HJ Sample 1 shows stable behavior while sample 2 eventually 
stops after running with lower energy for some time. Sample 3 shows larger fluctuation 
with higher energy; the sudden change of E(t) in sample 3 results from the change in the 
configuration of the bottom layer. This large fluctuation should be averaged out if we could 
simulate large enough system for long enough time, but the amount of computation is too 
large especially for the system with stiff particles. 

In order to make meaningful comparison out of these largely fluctuating data with a 
variety of behaviors, we select simulation sequences that come from similar flowing behaviors 
in the following way. First, we define the steady part of the time sequence in each of the 
samples as the part where the width of the energy fluctuations is smaller than 0.45 over the 
time period longer than 500. Second, we exclude the data whose averaged energy is out of 
the range [0.4,0.8] |22|. 

Many of the excluded data by this criterion show quite different flowing behaviors. We use 
only the data selected from this criterion to calculate time averages of physical quantities. 

In the snapshot of the frictional flow, Fig. most of the particles seem to be in contact 
with each other and form a layer structure. The y-dependence of the average density and 
the flow speed are shown in Figs. M (a) and (b), respectively. We can see that the relative 
motion between layers is largest at the bottom and very small in the bulk. The stiffness 
dependence of the average kinetic energy E is shown in Fig. El (a); the data are scattered 
due to the non-uniqueness of the steady state. 

Nevertheless, the average collision rates show systematic dependence as shown in Fig. 
ITUl (b). Here, the definition of N c (t c ) is the same as that for the collisional flow, namely, 
it is the average collision rate (contact time fraction) between particles per particle in the 
system. On the other hand, N w (t w ) is defined differently; it is the collision rate (the contact 
time fraction) between the particles and the slope per particle in the bottom layer, because 
other particles never touch the slope. 
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In Fig. El (b), we can see that N c (filled circles) and N w (open circles) increase very 
rapidly as k n becomes larger: they diverge as a power of k n . This is quite different from 
the behavior in the collisional flow in which they are almost constant. Furthermore, this 
increase is faster than the logarithmic divergence found in the single particle case (see Fig. 
13(b)). 

The contact time fractions, t c (filled circles) and t w (open circles), decrease as shown in 
Fig. El (c). The main reason why t c and t w decrease is that r c and r w decrease faster than 
N c and N w , namely, the contact time fractions estimated from the duration of a binary 
collision, N c t c and N w t w , continue to decrease. Actually, the decrease of the contact time 
fraction and the increase of the collision rate are natural because a longer multiple contact 
breaks up into shorter binary contacts as the particles become stiffer. 

These contact time fractions, t c and t w , however, do not converge to N c t c and N w t w , 
respectively. As shown in Fig. El (d), t c — N c r c (filled circles) and t w — N w r w (open 
circles) are very large as compared to those in the collisional flow, even in the stiffest region. 
The comparison of this with the rapid convergence in the collisional flow regime (Fig. 
(d)) indicates that there remains finite multiple contact time in the hard sphere limit and 
the interaction in the frictional flow can never be considered as many, or even infinite, 
instantaneous binary collisions. The particles experience the lasting multiple contact even 
in the hard sphere limit. 

IV. SUMMARY AND DISCUSSION 

The inelastic hard sphere limit of granular flow has been investigated numerically in the 
steady states of (i) a single particle rolling down the slope, (ii) the dilute collisional flow, and 
(iii) the dense frictional flow. In (i), it has been found that the "inelastic collapse" between 
the particle and the slope occurs in the hard sphere limit due to gravity, and the contact 
time fraction between the particle and the slope remains finite. In (ii), the collision rates 
N c and N w are almost constant when particles are stiff enough. The contact time fraction 
between particles t c approaches zero as k n increase in the same manner with the duration 
of contact for binary collision, r c . This means that the interaction between particles in the 
hard sphere limit can be expressed by binary collisions in the inelastic hard sphere model. 
On the other hand, the decrease in t w is slightly slower than l/yA^ in the harder region, 
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which can be a sign of the inelastic collapse between a particle and the slope. 

In the case of the frictional flow (iii), the situation is not simple: Although the contact 
time fractions t c and t w decrease upon increasing k n , the collision rates N c and N w increase 
as a power of k n , which is faster than the logarithmic divergence found in the single particle 
case (i). The origin of this power divergence of collision rate does not seem to be simple 
because it is a property of the steady state, not a particular dynamical trajectory of the 
system. 

The multiple contact time fraction may be estimated by t c — N c r c and t w — N w r w . They 
were found to be quite large compared to those for the collisional flow, and seems to remain 
finite even in the infinite stiffness limit: This suggests that the interaction in the frictional 
flow can never be considered as infinite number of binary collisions. Even in the hard sphere 
limit, particles experience the lasting multiple contact. 

The non-negligible fraction of multiple contact in the hard sphere limit implies the exis- 
tence of the network of contacting grains even though they are flowing. The models of dense 
flows should consider the effect of this lasting contacts. 

Here we have investigated only the two cases, i.e., the collisional flow and the frictional 
flow. The system should undergo the transition between the two flows if we change con- 
tinuously the parameters such as the inclination angle, the roughness of the slope, or the 
density of particles. It is interesting to investigate how the transition occurs by looking at 
the quantities measured in this paper, because their stiffness dependences are qualitatively 
different in the two flows. 

Acknowledgments 

This research was partially supported by Hosokawa powder technology foundation, the 
Japan Society for the Promotion of Science, and Grant-in- Aid from the Ministry of Educa- 
tion, Science, Sports and Culture of Japan. 



[1] S. B. Savage, Adv. Appl. Mech. 24, 289 (1984); C. Ancey, Phys. Rev. E 65, 011304 (2002). 
[2] J. T. Jenkins and S. B. Savage, J. Fluid Mech. 130, 187 (1983); C. S. Campbel, Annu. 
Rev. Fluid Mech. 22, 57 (1990); I. Goldhirsch, in Modeling in Applied Sciences: A Kinetic 



13 



Theory Approach, ed. N. Bellomo and M. Pulvirenti (Birkhauser, Boston, 2000); N. Mitarai, 
H. Hayakawa, and H. Nakanishi, Phys. Rev. Lett. 88, 174301 (2002). 
[3] P. Mills, D. Loggia, and M. Tixier, Europhys. Lett. 45, 733 (1999); B. Andreotti and S. 

Douady, Phys. Rev. E 63, 031305 (2001). 
[4] J. Duran, Sands, Powders, and grains: an introduction to the physics of granular materials 

(Springer, New York, 1997). 
[5] M. Isobe, Int. J. Mod. Phys. CIO, 1281 (1999). 
[6] M. Y. Louge, Phys. Fluids 6, 2253 (1994). 
[7] L. P. Kadanoff, Rev. Mod. Phys. 71, 435 (1999). 
[8] B. Berne and R. Mazighi, J. Phys. A: Math. Gen. 23, 5745 (1990). 
[9] S. McNamara and W. R. Young, Phys. Fluids A 4, 496 (1992). 
[10] S. McNamara and W. R. Young, Phys. Rev. E 50, R28 (1994); 53, 5089 (1996); N. Schorghofer 

and T. Zhou, ibid. 54, 5511 (1996); M. Alam and C. M. Hrenya, ibid. 63, 061308 (2001). 
[11] P. A. Cundall and O. D. L. Strack, Geotechnique 29, 47 (1979); A. Shimosaka, in Funtai 
Simulation Nyumon, edited by The Society of Powder Technology, Japan (Sangyo Tosho, 
Tokyo, 1998) p. 34, Sec. 3.2. 
[12] T. Pdschel, J. Phys. II (France) 3, 27 (1993). 
[13] N. Mitarai and H. Nakanishi, J. Phys. Soc. Jpn. 70, 2809 (2001). 

[14] D. Ertas et ai, Europhys. Lett. 56, 214 (2001); L. E. Silbert, G. S. Grest, S. J. Plimpton, and 
D. Levine, Phys. Fluids 14, 2637 (2002); L. E. Silbert et al, Phys. Rev. E 64, 051302 (2002). 
[15] S. Yuu, T. Abe, T. Saitoh, and T. Umekage, Adv. Powder Technol. 6, 259 (1995). 
[16] S. Luding, Phys. Rev. E 52, 4442 (1995). 

[17] J. Schafer, S. Dippel, and D. E. Wolf, J. Phys. I (France) 6, 5 (1996). 

[18] P. Henrici, Discrete Variable Methods in Ordinary Differential Equations (Wiley, New York, 
1962). 

[19] S. Dippel, G. G. Batrouni, and D. E. Wolf, Phys. Rev. E 54, 6845 (1996). 

[20] The stiffness dependence of the motion of particles discharging from a hopper using discrete 
element method has been studied by Yuu et al r|. They performed the simulations with three 
different values of stiffness and found the different behaviors. The comparison with the present 
results is not simple because of the different set up, but if we normalize their parameter of 
normal stiffness with the particle mass, the diameter, and the acceleration of gravity, then 



14 



their values extend over the range where E depends on k n in our simulations. 
When we fix one of the particles on the slope and employ the same force law as Eq. (J7J), then 
we have the duration of contact for the normal collision without gravity t w = ir/ \Jk n — rj^/A 
for identical particles. The force law for the flat boundary that replace 2M in Eq. (J7J) by the 
mass of the colliding particle m gives the same duration of contact with the above one for the 
normal collision. 

As for the softest case, l/k n = 2 7 x 10~ 5 , we had to include the data whose energy goes below 
the lower limit because of the general tendency that the energy becomes small for the steady 
part of the very soft region. For l/k n < 10~ 5 , however, we did not see any general trend in 
the average energy. 
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FIG. 1: A snapshot of a single ball rolling down a rough slope. 
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FIG. 2: Time evolution of v y with k n = 2 1 x 10 5 (solid line) and 2 17 x 10 5 (dashed line). 
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FIG. 3: The stiffness dependence of (a) the time averaged kinetic energy of the particle E, (b) the 
collision rate between the slope and the particle N w , and (c) the contact time fraction between the 
slope and the particle t w . 




FIG. 4: A snapshot of the dilute collisional flow. 
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FIG. 5: The stiffness dependences of (a) the averaged kinetic energy per one particle E, (b) the 
averaged collision rates between particles iV c (filled circles) and between particles and the floor 
N w (open circles), (c) the averaged contact time fractions between particles t c (filled circles) and 
between particles and the floor t w (open circles), (d) the estimated multiple contact time fractions, 
t c — N c t c (filled circles) and t w — N w t w (open circles). The solid and the dashed lines in (c) are 
proportional to 1/ \fh^. 
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FIG. 6: The number density profile (a) and the velocity profile (b) of the collisional flow. 
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FIG. 7: A snapshot of the dense fractional flow. 
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FIG. 8: Time evolution of the kinetic energy per particle of three samples (sample 1: red solid 
line, sample 2: green dashed line, sample 3: blue dotted line), (b) is the magnification of (a). 
Sample 1 with k n = 2 _1 x 10 5 shows steady behavior within the threshold. Sample 2 and 3 are 
with k n = 2~ 7 x 10 5 . Sample 2 shows steady behavior for a while but finally stops. E(t) of sample 
3 shoots up suddenly at t ~ 1300 when one of the particles in the bottom layer runs on other 
particles. 
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FIG. 9: The density (a) and velocity (b) profiles of the frictional flow. 
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FIG. 10: The stiffness dependences of (a) the averaged kinetic energy E, (b) the averaged collision 
rates between particles N c (filled circles) and between particles and the floor N w (open circles), (c) 
the averaged contact time fractions between particles t c (filled circles) and between particles and 
the floor t w (open circles), and (d) the estimated multiple contact time fractions, t c — N c t c (filled 
circles) and t w — N w t w (open circles). In (b) the solid line proportional to k^ A is shown for the 
guide of the eyes. 
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